#ifndef _GAME_SOLVER_HH_
#define _GAME_SOLVER_HH_
#include "matrix.h"
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <ostream>


std::ostream & dump(std::ostream & out, Matrix<double> & m)
{
  signed nrows = m.rows();
  signed ncols = m.columns();
  
    
  for ( int row = 0 ; row < nrows ; row++ ) {
    out << row << ':';
    for ( int col = 0 ; col < ncols ; col++ ) {
      out << " ";
      out.width(2);
      out << m(row,col);
    }
    out << std::endl;
  }
  out << std::endl;

  return out;
}



double sigmoid(double x)
{
  return 1.0/(1.0+exp(-x));
}


///
/// find optimal randomized strategy for both players
/// 
double solve(Matrix<double> & payoff, 
	     Matrix<double> & p, Matrix<double> & q, unsigned nmax)
{
  double lambda = 0.1;

  assert(payoff.minsize()>0);

  // makes the payoff matrix positive
  double minval = std::numeric_limits<double>::max();
  for ( int row = 0 ; row < payoff.rows() ; row++ )
    for ( int col = 0 ; col < payoff.columns() ; col++ ) 
      if(payoff(row,col)<minval) minval = payoff(row,col);

  for ( int row = 0 ; row < payoff.rows() ; row++ )
    for ( int col = 0 ; col < payoff.columns() ; col++ ) 
      payoff(row,col) -= minval;
   


  p.resize(1,payoff.rows());
  for(signed j=0;j!=p.columns();++j)
    p(0,j) = 1.0/static_cast<double>(p.columns());
  
  q.resize(payoff.columns(),1);
  for(signed i=0;i!=q.rows();++i)
    q(i,0) = 1.0/static_cast<double>(q.rows());


  double prev_reward = 42;
  double reward = p.product(payoff).product(q)(0,0);
  unsigned n = 0;
  for(n=0; n!=nmax and fabs(prev_reward - reward)>std::numeric_limits<double>::epsilon(); ++n)
    {
      Matrix<double> AliceGain = payoff.product(q);
      Matrix<double> BobGain = p.product(payoff);

      // reinforcement for Alice
      {
	for(signed i=0;i!=AliceGain.rows(); ++i)
	  p(0,i) *= exp(lambda*AliceGain(i,0));

	double Z=0;
	for(signed j=0;j!=p.columns(); ++j)
	  Z += fabs(p(0,j));
	for(signed j=0;j!=p.columns(); ++j)
	  p(0,j) /= Z;
      }


      // reinforcement for Bob
      if(0)
      {

	for(signed j=0;j<BobGain.columns(); ++j)
	  q(j,0) *= exp(-lambda*BobGain(0,j));

	double Z=0;
	for(signed i=0;i!=q.rows(); ++i)
	  Z += fabs(q(i,0));
	for(signed i=0;i!=q.rows(); ++i)
	  q(i,0) /= Z;
      }


	
      prev_reward = reward;
      reward = p.product(payoff).product(q)(0,0);
    }
 
  std::cerr << "Iterations: " <<  n << " steps."
	    <<std::endl;



  return reward;
 }


#endif /* _GAME_SOLVER_HH_ */

